Van der Waals (vdW) interaction typically refers to the London dispersion force, which is often qualitatively described as the interaction between an instantaneous dipole of a molecule and its induced dipole on another molecule. Sometimes Keesom interaction (i.e., thermally averaged interaction between two permanent dipoles) and Debye interaction (i.e., thermally averaged interaction between a permanent dipole and its induced dipole) are considered as generalized vdW interactions. All of these interaction energies follow a 1/r⁶ asymptotic behavior and the quantitative description of the vdW interaction requires quantum mechanics, which is different from classical electrostatic interactions. The London dispersion force is originated from pure non-local correlation of electrons, which is not captured by local or semi-local exchange–correlation functionals used in density functional theory (DFT).^[1,2] As a result, the DFT calculations using these functionals do not in principle contain the vdW interaction.^[3,4] Because DFT is the dominant approach for first-principles study on large systems, it has been an important topic to include the vdW interaction into DFT calculations.

Various methods have been proposed over the last two decades. Developing non-local vdW functionals has been pioneered by the Chamlers and Rutgers groups. These researchers developed a series of vdW-DF functionals.^[5–8] Later, Vydrov and van Voorthis proposed VV09 and VV10 functionals.^[9,10] Sabatini et al. presented a revision based on VV10 (rVV10), which allows the nonlocal correlation energy and its derivation to be evaluated in a plane wave framework while maintaining the precision of VV10.^[11] These vdW functionals are to be used together with semi-local functionals, for which significant progress has also been made over the years, such as TPSS,^[12] revTPSS,^[13] and M06L.^[14] In particular, the strongly constrained and appropriately normed (SCAN) meta-GGA functional obeys all 17 known exact constraints.^[15,16] Recently, the SCAN+rVV10 functional was extensively benchmarked and demonstrated high accuracy on layered compounds, where SCAN captures the short and intermediate-range interactions and rVV10 describes the long-range vdW interaction.^[17] Taking another philosophy, Grimme has developed pair potentials to be used together with the semi-local or hybrid functionals, dubbed DFT+D approach, which also demonstrates high accuracy and has received great attention in recent years.^[18–21] There are also other approaches in addition to the two approaches above, including dispersion-corrected atom-centered potentials (DCACP) and local atomic potentials (LAP).^[22–24]

The methods mentioned above are usually benchmarked by quantum chemistry methods, such as MP2 and CCSD(T), e.g., using the S22 database.^[25,26] The benchmark databases are usually in equilibrium structures and may not be suited for benchmarking the accuracy of DFT methods when sampling potential energy surfaces. For different sizes of systems, the methods or basis sets used in the benchmark databases are sometimes different, yielding inconsistency in the accuracy across the databases. It is desirable to develop consistent benchmark systems and systematically test the accuracy of the DFT methods, not only for the equilibrium configurations, but also for the overall potential energy surfaces.

In this paper, we present an accurate vdW benchmark system for potential energy surfaces, which were obtained with accuracy at the CCSD(T) level and complete basis set limit. Using this benchmark, we examined Grimme’s PBE+D3 method and the SCAN+rVV10 method, which were newly developed in recent years. The two methods represent the two most widely used treatments including the vdW interaction into DFT calculations, particularly for calculations on periodic systems using plane-waves as the basis set.^[27–36] In this paper, we focus on the plane-wave method, which is usually able to handle large systems and free of the basis set superposition error associated with local basis sets.^[37] In the next section, we will present the design of the benchmark system, which will be followed by comparison with the DFT methods, discussion on the results, and conclusion.

Our benchmark system is based on two molecules: glycine (an amino acid), denoted as Gly hereafter, and uracil (an RNA base), denoted as U hereafter. We considered three dimers formed by the two molecules, i.e., Gly-Gly, Gly-U, and U-U. For each dimer, we considered two potential energy surfaces (PESs). As illustrated in Fig. 1, one PES is with respect to the distance (d) and the other to the relative angle (θ ) between the two monomers. We thus have totally six PESs denoted by Gly-Gly-d, Gly-Gly-θ, Gly-U-d, Gly-U-θ, U-U-d, and U-U-θ, respectively. From Fig. 1, it can be seen that the three pairs represent three typical weak interaction systems. For the Gly-Gly pair, the London dispersion contributes more to the interaction energy than the other two pairs. The Gly-U pair is purposely arranged in a hydrogen-bonding configuration, which gives rise to stronger interaction than the other two pairs. The U-U pair represents a π–π stacking system, which usually yields stronger interaction than the pure London dispersion.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. (a) Gly-Gly dimer. (b) Gly-U dimer. (c) U-U dimer. Two small dots represent the centers of mass of the two monomers, which are connected by a dashed line. The arrows indicate the rotational direction.

The equilibrium distances for the d-PESs were taken to be the lowest energy points on the corresponding θ -PESs. For each PES, we sampled 20 points. The distance was measured between the centers of mass of the two monomers. For Gly-Gly-d, the distances (in units of Å) were taken as 2.7, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.9, 4.1, 4.3, 4.6, 4.9, 5.2, 5.7, 6.2, and 6.7. For Gly-U-d, the distances were taken as 5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6.0, 6.1, 6.2, 6.4, 6.6, 6.8, 7.1, 7.4, 7.7, 8.2, 8.7, and 9.2. For U-U-d, the distances were taken as 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.4, 3.5, 3.7, 3.9, 4.1, 4.4, 4.7, 5.0, 5.5, 6.0, and 6.5. For the θ PESs, the angels were taken to be uniformly separated by 18°.

We first optimized the atomic structures of the Gly and U monomers by MP2 calculations with the cc-pVTZ Gaussian basis set. In the dimer calculations, the structures of the monomers were fixed. The atomic structures at all data points are given in the supporting information (SI). At each data point on the PESs, the MP2 calculations were carried out with aug-cc-pVTZ and aug-cc-pVQZ basis sets, respectively. Basis set superposition error was corrected by the counterpoise method.^[37] Then, the MP2 binding energy at the complete basis set (CBS) limit was obtained by extrapolation using the Halkier’s scheme.^[38–40] CCSD(T) calculations were carried out using the aug-cc-pVDZ basis set to estimate the so-called ΔCCSD(T) correction. Overall, the CBS total energy was obtained by

where

Here, A and B are the fitting parameters and X equals to 3 and 4 for TZ and QZ basis sets, respectively. The correlation energy in the MP2 calculations is defined by . The ΔCCSD(T) term is usually less sensitive to the basis set and the aug-cc-pVDZ basis set is considered to yield accurate results for this term. Our MP2 and CCSD(T) calculations were carried out by using the MOLPRO package.^[41]

Our DFT calculations were performed using the VASP program.^[42] Both PBE (with and without the D3 correction) and SCAN (with and without the rVV10 correction) exchange–correlation functionals were used. The ion cores were represented by PAW potentials.^[43] Plane-waves with a kinetic energy cutoff of 544 eV were used as the basis set. The dimers were calculated in a cubic supercell of 25 Å along each side. The monomers were also fixed at the MP2-optimized structures and calculated using the same supercell. Grimme’s D3 parameters were taken from the tabulated values, which have been integrated into VASP.^[20,21] The interaction energy was calculated by taking the total energy difference between the dimer and the two monomers.

The purpose of designing the benchmark is to include both equilibrium structures and those away from the equilibrium, which are calculated on the same footing, therefore having the same level of accuracy in terms of the computational method and basis set. We designed the benchmark to include the representative systems, an amino acid and an aromatic RNA base molecule, which produce three systems of London dispersion-dominate, hydrogen bonding-dominate, and π–π stacking-dominate interactions, respectively. It is expected that this benchmark is easily reused by the community. For this purpose, we have included the geometries for all the 120 data points in the SI. All the raw data used in calculating the final results through Eqs. (1)–(4) are also included in the SI. The final MP2+ΔCCSD(T) results at the CBS limit, hereafter referred to as Gly-U benchmark or simply benchmark, are displayed in Figs. 2 and 3.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. Comparison between PBE, PBE+D3, and the benchmark on (a) Gly-Gly-d, (b) Gly-U-d, (c) U-U-d, (d) Gly-Gly-θ, (e) Gly-U-θ, (f) U-U-θ . All dashed and solid lines are for guiding the eyes.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. Comparison between SCAN, SCAN+rVV10, and the benchmark on (a) Gly-Gly-d, (b) Gly-U-d, (c) U-U-d, (d) Gly-Gly-θ, (e) Gly-U-θ, (f) U-U-θ.

First, we evaluate the accuracy of the PBE+D3 method. Figure 2 compares the PBE+D3 results with our benchmark on the six PESs. To illustrate the effect of the D3 correction, we also show the results from pure PBE. For the Gly-Gly and U-U pairs, the PBE method significantly underestimates the interaction energy, as is well known. From the results on Gly-Gly-θ and U-U-θ, it can be seen that PBE yields nearly no binding between the monomers from 0° to 360°, even though the general trend from PBE appears to be similar to the benchmark. By adding the D3 correction, the interaction energy is significantly improved, not only on the trend of the binding, but also point-by-point comparison. Different from Gly-Gly and U-U, for the Gly-U pair, i.e., the hydrogen-bonding system, the PBE already yields reasonably accurate results. More quantitative comparison will be given below together with the SCAN+rVV10 results.

Next, we evaluate the accuracy of the SCAN+rVV10 method. The comparison with the benchmark results is shown in Fig. 3. Again, to illustrate the effect of the rVV10 correction we show the SCAN results in Fig. 3. Compared with the PBE results in Fig. 2, the SCAN calculation already yields significantly improved results for all the three pairs. For the Gly-U pair, the SCAN method performs nearly equally well to the SCAN+rVV10 method, suggesting that for a hydrogen-bonding system, the SCAN method itself is sufficient to describe the interaction. For the pure London dispersion Gly-Gly pair and the π-π stacking U-U pair, the addition of rVV10 does further improve the interaction energy, yielding excellent agreement with the benchmark.

Figure 4 compares the deviation (ΔE) of the PBE+D3, SCAN, and SCAN+rVV10 results from the benchmark. Also shown in Fig. 4 is the mean average deviation (MAD) of each method. In terms of MAD, the SCAN+rVV10 method performs the best for the Gly-Gly and U-U pairs. In most cases, the SCAN+rVV10 method yields an error of 0.2 kcal/mol or smaller except for the case of Gly-Gly-d, where the SCAN+rVV10 method tends to overestimate the interaction (i.e., yielding less repulsive interaction if the interaction energy is positive or more attractive interaction if the interaction energy is negative) when d is smaller than the equilibrium distance. In comparison, the PBE+D3 method also performs well for the Gly-Gly and U-U pairs except for the case of U-U-d, where the deviation is large when d is smaller than the equilibrium distance. For the hydrogen-bonding Gly-U pair, the PBE+D3 method yields the smallest MAD, outperforming the SCAN and SCAN+rVV10 methods. The advantage of the PBE+D3 method for the Gly-U case is mainly in the region when d is smaller than the equilibrium distance, where the other two methods overestimate the interaction more significantly.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. Comparison of deviations of the PBE+D3, SCAN, SCAN+rVV10 interaction energies from the benchmark on (a) Gly-Gly-d, (b) Gly-U-d, (c) U-U-d, (d) Gly-Gly-θ, (e) Gly-U-θ, (f) U-U-θ . The arrows in (a)–(c) mark the equilibrium distances of the benchmark. The mean average deviation (in units of kcal/mol) of each method from the benchmark is also shown in the figure.

In general, if one considers only the region when d is around or greater than the equilibrium distance, both PBE+D3 and SCAN+rVV10 yield deviations well below the so-called chemical accuracy, i.e., 1 kcal/mol. However, when d is smaller than the equilibrium distance, neither of these methods yield as good accuracy uniformly for all the London dispersion, hydrogen-bonding, and π–π stacking systems. Such regions might be important for studying systems under high pressure or drastic conditions, such as explosive reactions.

It is worth to mention that in the curve for Gly-Gly-θ, there is a detailed structure in the region from about 50° to 130°, for which the SCAN and SCAN+rVV10 methods can reproduce the bump at 90° and the dip at 126°, while the PBE+D3 method fails to reproduce these features. Actually, the PBE calculation without D3 already fails to capture these subtle features. The D3 method is not able to correct the failure of PBE.

Finally, as the local density approximation (LDA) is commonly used in studying vdW systems,^[44] we also evaluated this method based on the Gly-U benchmark. Figure 5 shows the comparison results. In general, LDA overestimates the interaction as it does for most bulk materials. This is particularly true for pure London dispersion interaction, as can be seen for the Gly-Gly pair. Interestingly, the LDA curve for Gly-Gly-d decays much faster than the benchmark as d increases, suggesting that the vdW interaction in LDA is of wrong origin and acts only in a short range. For the Gly-U pair, LDA also overestimates the interaction in the repulsive region as well as in the region near the equilibrium. Only in the region when the two monomers are farther separated, it yields similar results to the benchmark. While LDA works reasonably good for Gly-U-θ, the agreement with the benchmark is still visibly worse than PBE+D3 and SCAN+rVV10. Similar to PBE and PBE+D3, LDA also fails to reproduce the subtle feature in Gly-U-θ from 50° to 130°.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. Comparison between LDA results and the benchmark on (a) Gly-Gly-d, (b) Gly-U-d, (c) U-U-d, (d) Gly-Gly-θ, (e) Gly-U-θ, (f) U-U-θ.

It is well known that LDA can accidentally produce a good inter-layer distance for graphite, which can be considered as a π–π stacking system. For the U-U pair, LDA can indeed generate good interaction energy near the equilibrium. But, when d is larger than equilibrium, the interaction energy curve decays more quickly than the benchmark, similar to the case of Gly-Gly-d. It is interesting that LDA does not consistently overestimate the vdW interaction. For U-U-d, as well as for Gly-Gly-d, LDA could also underestimate the vdW interaction when d is large. For U-U-θ, while consistently underestimating the interaction, LDA appears to describe the π–π stacking better than the pure London-dispersion in Gly-Gly-θ.

In summary, we designed a benchmark system for the vdW interaction with an aim to test the whole potential energy surfaces instead of individual configurations only. Three systems, namely, Gly-Gly, Gly-U, and U-U, represent London dispersion, hydrogen bonding, and π–π stacking interactions, respectively. The benchmark was obtained with the MP2+ΔCCSD(T) method at the CBS limit. With this benchmark, we examined the PBE+D3 and SCAN+rVV10 methods. Both methods can yield accuracy well below 1 kcal/mol, except for the cases when the distance between two monomers is below the equilibrium distance. In such a case, neither of these methods is uniformly accurate for the London dispersion, hydrogen-bonding, and π–π stacking systems. In addition, it was found that the SCAN functional can reproduce some subtle features in a potential energy curve, while both PBE and LDA fail. When PBE fails to reproduce such features, the D3 method is not able to correct this failure.